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Abstract. We describe a new approach to gravitational instability in 
large-scale structure, where the equations of motion are written and 
solved as in field theory in terms of Feynman diagrams. The basic ob- 
jects of interest are the propagator (which propagates solutions forward in 
time), the vertex (which describes non- linear interactions between waves) 
and a source with prescribed statistics which describes the effect of initial 
conditions. We show that loop corrections renormalize these quantities, 
and discuss applications of this formalism to a better understanding of 
gravitational instability and to improving non-linear perturbation theory 
in the transition to the non- linear regime. We also consider the role of vor- 
ticity creation due to shell-crossing and show using N-body simulations 
that at small (virialized) scales the velocity field reaches equipartition, 
i.e. the vorticity power spectrum is about twice the divergence power 
spectrum. 



1. Standard Formulation of Gravitational Instability 

Assuming the initial velocity field is irrotational, gravitational instability can be 
described completely in terms of the density field and the velocity divergence, 
6 = V • v. Defining the conformal time r = J dt/a and the conformal expansion 
rate H = dlna/dr, the equations of motion in Fourier space become 



Equations (|T]) and (||) are valid in an arbitrary homogeneous and isotropic uni- 
verse, which evolves according to the Friedmann equations: 




where [<5d] = o~D(k — ^12), fe is a comoving wave number, and 




(3) 
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(n-l)H 2 (T) = k-^a 2 (r), 



(5) 



where A is the cosmological constant, the spatial curvature constant k = —1,0, 1 
for Q tot < 1, £l tot = 1, and O to t > 1, respectively, and £7 to t = ^ + ^A, with 
J7a = Aa 2 /(3W 2 ). For Q = 1, the perturbative growing mode solutions are given 
by 



S(k) = £ a n (r)6 n (k), 

n=l 

oo 

0(k) = W(r) $> n (r)0 n (k). 



(6) 
(7) 



n=l 



Modelling the matter as pressureless non-relativistic 'dust', an appropriate de- 
scription for cold dark matter before shell crossing, the fluid equations of motion 
determine 8 n {k) and n (k) in terms of the linear fluctuations, 

5 n (k) = Jd 3 q 1 ...Jd i q n [5 D ] n F r [ s \q 1 ,...,q n )6 1 (q 1 )...8 1 (q n ), (8) 
On(k) = ~Jd 3 q l ...Jd 3 q n [8 D ] n G l n s \q 1 ,...,q n )5 l (q 1 )...5 1 (q n ), (9) 

where [<5d]« = Sr>(k — q^ — . . . — q n ). The functions Fn and g4 are constructed 
from the mode coupling functions a(k, fei) and /?(fci, £2) by a recursive procedure 
(Goroff et al. 1986), 



Fn{q\, ■■■,q n ) 



G n (qi, ■ ■ ■ ,q n ) 



n-1 

E 



Gm(qi, ■ ■ ■ > q Tl 



f x (2n + 3)(n-l 
+2/?(fc 1 ,fc 2 )G 
G m (qi, ■ ■ ■ ,q s 



(2n + l)a(fe,fci)F„_ m ( 

"m+l > 



E 



3a(fc,fei)i ? n 

-m(9m+b • • • j 9n) 



'j (2n + 3)(n-l) 
+2n/3(fci, k2)G n - m (q m+ i, . . . ,q n ) 



,q n ) 

(10) 

(ii) 



(where fei = q t + . . .+q m , k 2 = q m+ i + . . • +<?„, k = ki+k 2 , and F\ = G\ = 1). In 
Eqs.dli), Fl s) and G$ are the symmetrized version of F n and G n , respectively. 



From these perturbative solutions a number of important results have been 
derived in the literature, most of them regarding the tree-level (leading-order) 
behavior of correlation functions, e.g. Fry (1984), Goroff et al. (1986), Bernardeau 
(1992,1994). Loop calculations have been attempted only in some particular 
cases (Scoccimarro & Frieman 1996; Scoccimarro 1997, Scoccimarro et al. 1998); 
although in the spherical collapse approximation a number of useful results have 
been obtained (Fosalba & Gaztahaga 1998; Gaztahaga & Fosalba 1998). 
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2. Field Theory Approach 



2.1. Integral Form of the Equations of Motion 

The equations of motion can be rewritten in a more symmetric form by defining 
a two-component "vector" ^ a (k,z), where a = 1, 2, z = lna, and: 

*a(*s, z) = \S(k, z), -9{k, z)/n\ , (12) 

which for 17 = 1 leads to the following equations of motion (we henceforth use 
the convention that repeated Fourier arguments are integrated over) 

d z ^ a (k,z) + Q ab ^ b (k,z) 

where 

and j a b c is a matrix given by 

7i2i(fc,fci,fc2) = <5 D (fc - fci2) a(k,ki), (15) 
7222(^,^1,^2) = S D (k - k 12 ) /3(fci, fc 2 ), (16) 

and zero otherwise. The somewhat complicated expressions for the PT kernels 
recursion relations in the previous section can be easily derived in this formalism. 
The perturbative solutions read 

00 

* a (fc,z) = £e n2 Vi n) (fc), (17) 

71=1 

which leads to 

n-l 

(n5 ab + n ab ) 4 n \ k )=7abc(k,ki,k 2 ) £ Vi n " m) (fci) ^ m) (fe 2 ). (18) 

m=l 

Now, let (?~ b {n) = n8 ab + £l ab , then we have: 

n-l 

^W(fc) = a ab (n) 76c <i(fe,fel,fe2) E ^" m) (fci) ri m) (fc 2 ), (19) 

m=l 

where 

°"afc( n ) = 77^ T^TT TT 

(2n + 3)(n — 1) 

Equation ( |i"9| ) is the equivalent to the Goroff et al. (1986) recursion relations, 
Eqs. ( |iTj[Ji"l|) , for the n th order Fourier amplitude solutions V'a^(fc). It turns out 
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Jabc{k,k 1 ,k 2 ) ^b{ki,z) y c {k 2 ,z), (13) 



-1 
-3/2 1/2 



(14) 



2n + 1 2 
3 2n 



(20) 



to be convenient to write down the equation of motion Eq. (13) in integral form. 
Laplace transformation in the variable z leads to: 



®b(k,u) = (t>a(.k)+j abc (k,ki,k 2 ) I ^ b (k 1 ,u 1 )^ c (k 2 ,u; -Ui), (21 



where (j) a {k) denote the initial conditions, that is ^ a (k,z = 0) = (f> a {k). Multi- 
plying by the matrix o~ ab , and performing the inversion of the Laplace transform 
we finally get 



* (fe, z) = g a b(z) (f)b(k) + / dz' g ab (z - z) j b cd(k, k l: k 2 ) ^ c (ki,z')^d(k 2 , z'), 

Jo 

(22) 

where the linear propagator g a b(z) is defined as (c > 1 to pick out the standard 
retarded propagator, Scoccimarro 1998) 



g a b(z) 



c+ioo duJ 
:< 

c— zoo 27TZ 



e 2 


3 2" 


e -3z/2 


-2 2 " 


~5~ 


3 2 


5 


3 -3 



(23) 



for z > 0, whereas g a b( z ) = for z < due to causality, g ab (z) — » S a b as 
z — > + . The first term in Eq. ( |2"3| ) represents the propagation of linear growing 
mode solutions, where the second corresponds to the decaying modes propaga- 
tion. If we assume that the initial conditions are set in the growing mode, then 
4>a(k) = 5i(fe)(l, 1) and the second term in Eq. (23) does not contribute upon 
contraction with (j>b(k)- Consistently with this, we can neglect subdominant time 
dependences in the non- linear term in Eq. (22), which amounts to setting the 
initial conditions at z = —00. Then, the equations of motion in integral form 
reduce to: 



^ a (k,z) = e z M k )+ I dz' g a b(z-z')~fbcd(k,k 1 ,k 2 )y c (k ll z')y d (k2,z'). (24) 

^—00 

As it stands, this integral equation can be thought as an equation for x I / a (fc, z) in 
the presence of an "external source" (j) a (k) with prescribed statistics. In partic- 
ular, if we assume that the initial conditions are Gaussian; then the statistical 
properties of 4> a (k) are completely characterized by its two-point correlator 

( Mk) Mk') )=S D (k + k') u ab P(k), (25) 

where P(k) denotes the initial power spectrum of density fluctuations and u a b = 
1 for growing-mode initial conditions. From Eq. (24), it is easy to verify that 
the ansatz in Eq. ( |i~7| ) leads to the recursion relations in Eq. (|l~9|). 

Equation (22) has a simple interpretation. The first term corresponds to the 



linear propagation from the initial conditions, whereas the second term contains 
information on non- linear interactions (mode-mode coupling). This corresponds 
to all the interactions between waves that happen at time z' (with < z' < z) 
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characterized by 7^ and then propagated forward in time from z' to z by 
the propagator g ab {z — z'). We can also write down the general solution as a 
perturbation series, 



^(k,z) = J 5 D (k-k 1 ... n )^ n \k 1 ,...,k n ;z)S 1 (k 1 )...5 1 (k n ), (26) 
where the kernels satisfy the usual recursion relations 



n fZ 

^ n) (*)=E / dsg ah {z-s) lhcd {kMM)T^\s)T { r m \s). (27) 

m=l J ° 

Interactions modify the linear propagator, leading to propagator renormal- 
ization, so that the non-linear propagator defined by 



G ab (k,z) 5b (k - k!) = 



6^ a (k,z) 
Hb(k') 



reads 



00 „ 

Gab{k, z) = g ab {z) + A n (z) I cPq^ . . . d?q n P, 

n=l 



n n 



(2n+l) 



(28) 



(29) 



=(2n+l) 



where F£'*^> = ^ 2n+1) (k,q 1: -q x , ...,q n , -q n ), A n (z) = (2n - l)!!exp[(2n + 
l)z], and we defined (f) b = (ui,U2)5i(k). Similarly the vertex is renormalized by 
non-linear interactions as well, 



Fabc(ki,k 2 ,z) 6b (k - k 12 ) = G bd G C g 



6 2 * a (k,z) 



(30) 



and thus T abc = ^ abc + corrections. 

The calculation of correlation functions can be written down in terms of 
Feynman diagrams (Figs. We assign a solid line to each propagator, 

Eq. (p3|), a crossed circle represents the two-point correlator in the initial con- 
ditions, Eq. (25), and a solid circle represent the vertex, Eqs. (15-0|). In this 
representation, loop corrections can be divided into two general classes, those 
which renormalize the propagator, and those which renormalize the vertex. For 
example, in the calculation of the one-loop power spectrum there are two con- 
tributions, = {6163 + 6262)0 (where <5j denotes the i th order solution in 
PT); the (<5i<$3) c term corresponds to renormalizing the propagator (first one- 
loop term in Fig. |]), whereas the {6262 ) c term denotes the irreducible one-loop 
power spectrum (second one- loop term in Fig. |l]). By irreducible, we mean that 
this contribution cannot be separated into two connected diagrams by cutting 
one internal propagator line, unlike the ( 5\5$ ) c contribution. 

For the bispectrum, the 4 one-loop terms can be divided in a similar fashion. 
The ( 6461 ) c corresponds to vertex renormalization (first one-loop term in Fig. Q), 
the ( 636261 ) c correspond to power spectrum (second one-loop term in Fig. 0) and 
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propagator renormalization (third one-loop term in Fig. Q), and the (<5§ } c gives 
the irreducible one-loop bispectrum (last term in Fig. ||) . This formalism can also 
be extended to include non-Gaussian initial conditions, see Scoccimarro (1998) 
for a general discussion and the specific example of Zel'dovich approximation 
initial conditions, relevant to transients in N-body simulations. 

3. One-Loop Propagator and The Non-Linear Power Spectrum 

As an example, we calculate the one-loop propagator, G^j 

G { V(k,z)=g ab (z)+ex V (3z) J d 3 qP(q)-^-(k, q, -q; z), (31) 

If we take the k — > limit, we find that (keeping only the fastest growing 
term) 

G ab (&> z ) = 9ab{z) - o% exp(3z) 

where a 2 = J P(q)d 3 q/q 2 . Since the correction is negative, this tends to make 
the non-linear growth smaller than in linear theory, particularly for linearly 
decaying modes, which decay faster than in the linear case. The correction to 
g ab can be rewritten in terms of a correction to Q ab , using that 

(d z G ab ) Gfc = -Vac, (33) 

so that 

5n ab « k 2 a 2 v exp(7z/2) 

In order to see the role of decaying modes in the standard solutions of non- 
linear PT, let's consider the usual second-order PT kernel (e.g. relevant for the 
calculation of the skewness and bispectrum). In our notation, the second-order 
kernel can be written as 



9/50 61/1050 
3/25 61/1575 



(32) 



-28/375 28/375 
-4/75 4/75 



(34) 



J r 1 (2) (fci,fc 2 ) = / ds g u (z - s) -y bcd exp(2s) (1, l) c (1, l) d , (35) 
Jo 

where we assumed linear growing mode initial conditions. Since g\ b j bc d = 
5n7i2i + 5127222 1 arid using Eq. (f23|) we have for the fastest growing contri- 
bution to 2nd-order PT 



^"l (2) (fel,fe2) 



exp(2z) 



(^ + 4) a(fc ' fcl) + (^-4) /3(fcl ' fc2) 



r5 2 
exp(2z) y-a{k, fci) + -/3(fei, fe 2 ) 



(36) 
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It is crucial to note here that the 4/35 contributions come from linearly decaying 
modes; that is, after scattering, waves are not in linearly growing modes any- 
more, and this type of amplitude propagated into the present time contributes 
4/35 to the amplitudes of second-order PT kernels. That means that if we 
are using the kernels to calculate one-loop corrections, which are important at 
intermediate k's, one could use the approximation in Eq. (|34|) to improve the 
propagator in Eq. (|35|); in this case this corresponds to supressing the linearly 
decaying mode contribution to the propagator. As a result, the 2nd-order PT 
kernel at intermediate scales would look more like 



This means that a simple way of improving one-loop corrections in PT, is to use 
the kernels obtained by ignoring the linearly decaying modes contributions. This 
has the effect of incorporating higher-order loop corrections (those corresponding 
to propagator renormalization, although only approximately since we don't use 
the full one-loop propagator) in the usual formulation of PT. 

If we supress linearly decaying mode contributions for the third-order kernel 
and use this to calculate one-loop corrections to the power spectrum, we find 
the results in Figs||]|. The solid lines in Fig. ^ show the standard (top) and 
"improved" (bottom) calculations of the non-linear power spectrum, whereas 
the dashed line shows the fitting formula for the non-linear power spectrum. 
The three-remaining solid lines (which extend up to k ~ 3 h/Mpc) denote the 
measurement in N-body simulations of the density power spectrum, the velocity 
divergence power spectrum and the velocity vorticity power spectrum, as labeled. 
We see that the vorticity power spectrum is certainly negligible at large scales, 
and it does not become significant until scales of order k ~ 2 h/Mpc. At small 
scales, we find that the vorticity spectrum is roughly twice that of the divergence, 
as expected if the velocity field has equal power in all directions relative to k. 

Note that the "improved" calculation is somewhat smaller than the stan- 
dard one-loop calculation, as expected since the contribution from propagating 
linearly decaying modes has been suppressed. Overall the agreement with the 
N-body results is better. In Fig. |] we show the ratio of our predictions for dif- 
ferent models to the non-linear fitting formula, the horizontal dashed lines show 
the expected accuracy of the latter. We see that the "improved" calculations 
(solid) stay within the non-linear fitting formula accuracy up to k ~ 5 h/Mpc, 
whereas the standard one-loop calculation (dashed) overestimates the non-linear 
power spectrum at scales smaller than the non-linear scale, k ~ 0.3 h/Mpc. The 
improvement is thus quite significant, although we have included the effects of 
propagator renormalization in a crude way. 

Unlike the velocity divergence, which can be calculated in one-loop PT in 
analogous fashion to the density power spectrum, understanding the vorticity 
power spectrum is considerably more complicated because vorticity is generated 
by shell-crossing, an effect which is negelected in the formulation of PT (see e.g. 
Pichon & Bernardeau, 1999). However, we can understand approximately the 
scaling with redshift and scale from simple considerations. After shell-crossing, 
vorticity develops because what we see is the mass average of different streams, 
each with its own (irrotational) velocity field. Thus, vorticity can be thought as 
coming from the vorticity of the mass-weighted velocity field, i.e. 




(37) 
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w~f v (r) Vx l(l + S)v], (38) 

where f v {r) is the fraction of volume that undergoes shell-crossing at time r. We 
can then write the vorticity power spectrum (( w(k) ■ w(k') ) = P w (k)5i)(k + k')) 

P w (k) ~ f v {T) J ^^[P e (\k - q\)P s (q) - P x (\k - q\)P x (q)]d 3 q, (39) 

where Pe(k) is the velocity divergence power spectrum and P x (k) is the power 
spectrum of the density-velocity divergence cross-correlation. The simplest ap- 
proximation would be to use linear PT (although it is unlikely to be valid for 
each flow at the scales of interest); however, since Pg = P$ = P x in linear PT, 
this contribution vanishes. Thus, the leading-order contribution to P w (k) comes 
from one-loop PT, 



P w (k) ~ f v (r)k 2 J ^a*\k - g| 2 "+V ~ a 6 /^ 3ra+6 . (40) 



d 3 q 
"? ' 

Thus, we expect a strong time and scale dependence for the vorticity power spec- 
trum. The latter is in reasonable agreeement with Fig. ||, the time dependence 
is more difficult to test due to the unknown dependence coming from f 2 (r). 



4. Conclusions 

We described a new approach to gravitational instability in large-scale struc- 
ture, where the equations of motion are written and solved as in field theory 
in terms of Feynman diagrams. The basic objects of interest are the propaga- 
tor (which propagates solutions forward in time), the vertex (which describes 
non-linear interactions between waves) and a source with prescribed statistics 
which describes the effect of initial conditions. Loop corrections renormalize 
these quantities, in particular, decaying modes are supressed in the one-loop 
propagator compared to linear PT. We used this to construct the PT kernels 
and calculate "improved" loop corrections, these include effects beyond standard 
one-loop PT calculations, leading to better agreement with N-body simulations 
for the evolution of the power spectrum. We also consider the role of vortic- 
ity creation due to shell-crossing and show using N-body simulations that at 
small (virialized) scales the velocity field reaches equipartition, i.e. the vorticity 
power spectrum is about twice the divergence power spectrum. We also sketched 
a derivation of the time dependence and scaling of the vorticity power spectrum. 

Our calculations are only a first attempt to include the effects of propagator 
renormalization, more work is needed to confirm that the results presented here 
are indeed robust to a more careful treatment. We have also neglected vertex 
renormalization. On the other hand, it seems that this approach can lead to 
useful insights into the nature of non-linear corrections and perhaps give us a 
more accurate way to calculate clustering statistics in the transition to the non- 
linear regime. Our results on the vorticity from numerical simulations suggest 
that there is a significant range of scales until the assumption of irrotational fluid 
breaks down. We hope that by using these techniques we can finally answer the 
question: 'Why does PT work so well?" 
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Figure 1. Power spectrum diagrams up to one-loop. The first term 
denotes the linear contribution, the two remaining terms denote the 
one-loop correction. The factor enclosed by dashed lines denotes prop- 
agator renormalization. 




Figure 2. Bispectrum diagrams up to one-loop. The first terms de- 
notes the tree contribution, the four remaining terms the one-loop cor- 
rection. The first factor enclosed by dashed lines denotes vertex renor- 
malization, the second corresponds to the irreducible one-loop power 
spectrum, the third denotes propagator renormalization. The last term 
gives the irreducible one-loop bispectrum. 
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k [h Mpc-'] 

Figure 3. The power spectrum of the density, velocity divergence 
and vorticity as a function of scale. The two solid lines roughly parallel 
at high-k are the standard one-loop PT density power spectrum calcu- 
lation (top) and the new one-loop approach (bottom). The dashed line 
denotes the prediction of the fitting formula and the solid line close to 
it the actual measurement in the N-body simulation. The two other 
solid lines denote the power spectrum of the velocity divergence and 
vorticity, as labeled. 
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Figure 4. Ratio of predictions from one-loop PT (standard in dashed 
lines, new approach in solid lines) to the fitting formula for the non- 
linear power spectrum. The three different curves are for three different 
models, as labeled. 
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